// todo
// Rewrite C-like file handling to a C++ version.
//
// Some new things to improve...
// rewrite this file changing verts|materials|triangles variables to object3dv


// All of the functions assume that they're being called with sensful arguments.

#include "obj_loader.h"
#include "error.h"
#include "material.h"
#include "triangle.h"
#include "vertex.h"
#include "vect.h"
#include "object3d.h"

#include <iostream>
#include <map>
#include <vector>
#include <utility>
#include <cmath>
#include <list>
#include <algorithm>
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <assert.h>
using std::cout;
using std::cerr;
//using std::endl;

static const float INF = 1e10;
static const int TRIANGLES_CONSTRAINT = 10;

error load_gat(const char* path, std::vector<vertex> & verts, int & verts_count,
               std::vector<material>& materials, int& materials_count,
               std::vector<triangle>& triangles, int& trcount, std::pair<float,float>& xborders, 
               std::pair<float,float>& yborders, std::pair<float,float>& zborders)
{
    FILE* file;

    file = fopen(path, "r");
    if(file == NULL)
        return error(error::ERROR_couldnot_open_file,0);
    
    int count;
    char line[16];
    char* r;
    int i;
    int vcount;
    int act_line = 1;
    int normalcount;
    int materialcount;
    int tcount;

    // The first line should be just "vertices".
    r = fgets(line, 16, file);
    if(r!=line)
    {
        fclose(file);
        return error(error::ERROR_file_IO,act_line);
    }
    i=0;
    while(i<16)
    {
        if(line[i]=='\0')
            break;
        else
            i++;
    }
    if(i==15 && line[14]!='\n')
    {
        // We haven't read the whole line. 
        fclose(file);
        return error(error::ERROR_corrupted_data,act_line);
    }
    
    act_line++;
    count = fscanf(file, "%d", &vcount);
    if(count!=1 || vcount<=0)
    {
        fclose(file);
        return error(error::ERROR_corrupted_data,act_line);
    }
    
    act_line++;
    //todo
    verts.reserve(vcount);
    for(i=0; i<vcount; i++)
    {
        verts.push_back(vertex());
        count = fscanf(file, "%f %f %f", &verts[i].xyz[0], &verts[i].xyz[1], &verts[i].xyz[2]);
        if(count!=3)
        {
            fclose(file);
            return error(error::ERROR_corrupted_data,act_line);
        }
        if(verts[i].xyz[0] < xborders.first)
           xborders.first = verts[i].xyz[0];
        if(verts[i].xyz[0] > xborders.second)
           xborders.second = verts[i].xyz[0];

        if(verts[i].xyz[1] < yborders.first)
           yborders.first = verts[i].xyz[1];
        if(verts[i].xyz[1] > yborders.second)
           yborders.second = verts[i].xyz[1];

        if(verts[i].xyz[2] < zborders.first)
           zborders.first = verts[i].xyz[2];
        if(verts[i].xyz[2] > zborders.second)
           zborders.second = verts[i].xyz[2];
        act_line++;
    }

    
    act_line--;
    // Skip "white" remainings of a previous block.
    r = fgets(line, 16, file);
    if(r!=line)
    {
        fclose(file);
        return error(error::ERROR_file_IO,act_line);
    }
    act_line++;
    r = fgets(line, 16, file);
    if(r!=line)
    {
        fclose(file);
        return error(error::ERROR_file_IO,act_line);
    }

    act_line++;
    // The second block starts with "normal vectors".
    r = fgets(line, 16, file);
    if(r!=line)
    {
        fclose(file);
        return error(error::ERROR_file_IO,act_line);
    }
    i=0;
    while(i<16)
    {
        if(line[i]=='\0')
            break;
        else
            i++;
    }
    if(i==15 && line[14]!='\n')
    {
        // We haven't read the whole line. 
        fclose(file);
        return error(error::ERROR_corrupted_data,act_line);
    }


    act_line++;
    count = fscanf(file, "%d", &normalcount);
    if(count!=1 || normalcount<=0 || normalcount!=vcount)
    {
        fclose(file);
        return error(error::ERROR_corrupted_data,act_line);
    }
    
    act_line++;
    for(i=0; i<vcount; i++)
    {
        count = fscanf(file, "%f %f %f", &verts[i].normal[0], &verts[i].normal[1], &verts[i].normal[2]);
        if(count!=3)
        {
            fclose(file);
            return error(error::ERROR_corrupted_data,act_line);
        }
        act_line++;
    }


    act_line--;
    // Skip "white" remainings of a previous block.
    r = fgets(line, 16, file);
    if(r!=line)
    {
        fclose(file);
        return error(error::ERROR_file_IO,act_line);
    }
    act_line++;
    r = fgets(line, 16, file);
    if(r!=line)
    {
        fclose(file);
        return error(error::ERROR_file_IO,act_line);
    }


    act_line++;
    // The third block starts with "materials".
    r = fgets(line, 16, file);
    if(r!=line)
    {
        fclose(file);
        return error(error::ERROR_file_IO,act_line);
    }
    i=0;
    while(i<16)
    {
        if(line[i]=='\0')
            break;
        else
            i++;
    }
    if(i==15 && line[14]!='\n')
    {
        // We haven't read the whole line. 
        fclose(file);
        return error(error::ERROR_corrupted_data,act_line);
    }

    act_line++;
    count = fscanf(file, "%d", &materialcount);
    if(count!=1 || materialcount<0)
    {
        fclose(file);
        return error(error::ERROR_corrupted_data,act_line);
    }

    // Skip the remaining of the previous line.
    r = fgets(line, 16, file);
    if(r!=line)
    {
        fclose(file);
        return error(error::ERROR_file_IO,act_line);
    }
    
    
    if(materialcount>0)
    {
       //todo
        materials.reserve(materialcount);
        act_line++;
        for(i=0; i<materialcount; i++)
        {
            materials.push_back(material());
            count = fseek(file, sizeof("diffuse ")-1, SEEK_CUR);
            if(count!=0)
            {
                fclose(file);
                return error(error::ERROR_file_IO,act_line);
            }
            count = fscanf(file, "%f %f %f", &materials[i].diffuse[0], &materials[i].diffuse[1], &materials[i].diffuse[2]);
            if(count!=3)
            {
                fclose(file);
                return error(error::ERROR_corrupted_data,act_line);
            }

            // Move over the remaining of the last line.
            r = fgets(line, 16, file);
            if(r!=line)
            {
                fclose(file);
                return error(error::ERROR_file_IO,act_line);
            }

            act_line++;
            count = fseek(file, sizeof("specular ")-1, SEEK_CUR);
            if(count!=0)
            {
                fclose(file);
                return error(error::ERROR_file_IO,act_line);
            }
            count = fscanf(file, "%f %f %f", &materials[i].specular[0], &materials[i].specular[1], &materials[i].specular[2]);
            if(count!=3)
            {
                fclose(file);
                return error(error::ERROR_corrupted_data,act_line);
            }

            // Move over the remaining of the last line.
            r = fgets(line, 16, file);
            if(r!=line)
            {
                fclose(file);
                return error(error::ERROR_file_IO,act_line);
            }

            act_line++;
            r = fgets(line, 16, file);
            if(r!=line)
            {
                fclose(file);
                return error(error::ERROR_file_IO,act_line);
            }

            act_line++;
        }
        
        act_line--;
    }
    else
    {
        act_line++;
        r = fgets(line, 16, file);
        if(r!=line)
        {
            fclose(file);
            return error(error::ERROR_file_IO,act_line);
        }
    }

    act_line++;
    r = fgets(line, 16, file);
    if(r!=line)
    {
        fclose(file);
        return error(error::ERROR_file_IO,act_line);
    }
    i=0;
    while(i<16)
    {
        if(line[i]=='\0')
            break;
        else
            i++;
    }
    if(i==15 && line[14]!='\n')
    {
        // We haven't read the whole line. 
        fclose(file);
        return error(error::ERROR_corrupted_data,act_line);
    }
    
    act_line++;
    count = fscanf(file, "%d", &tcount);
    if(count!=1 || tcount<=0)
    {
        fclose(file);
        return error(error::ERROR_corrupted_data,act_line);
    }
    
    //todo
    triangles.reserve(tcount);
    act_line++;
    for(i=0; i<tcount; i++)
    {
        triangles.push_back(triangle());
        count = fscanf(file, "%d %d %d", &triangles[i].verts[0], &triangles[i].verts[1], &triangles[i].verts[2]);
        if(count!=3)
        {
            fclose(file);
            return error(error::ERROR_corrupted_data,act_line);
        }
        act_line++;
    }

    verts_count = vcount;
    materials_count = materialcount;
    trcount = tcount;
    
    fclose(file);
    return error(error::OK,0);
}


// I assume that object structure is correct, e.i.:
// verts.size() == verts_count
// materials.size() == materials_count
// triangles.size() == trcount
// verts_count >= 3 (at least one triangle exist)
// materials_count >= 0
// trcount >= 1
// there are no duplicate vertices neither triangles nor materials
int prepare_lod_base(std::vector<vertex>& verts, int verts_count,
               std::vector<material>& materials, int materials_count,
               std::vector<triangle>& triangles, int trcount,
               object3dv<tvertex>& lodbase,
               vect& mins, vect& maxs)
{
    int i;
    lodbase.vects.resize(verts_count);
    mins.x = INF;
    mins.y = INF;
    mins.z = INF;
    maxs.x = -INF;
    maxs.y = -INF;
    maxs.z = -INF;
    for(i=0; i<verts_count; i++)
    {
        lodbase.vects[i].xyz[0] = verts[i].xyz[0];
        if(verts[i].xyz[0] < mins.x)
           mins.x = verts[i].xyz[0];
        if(verts[i].xyz[0] > maxs.x)
           maxs.x = verts[i].xyz[0];

        lodbase.vects[i].xyz[1] = verts[i].xyz[1];
        if(verts[i].xyz[1] < mins.y)
           mins.y = verts[i].xyz[1];
        if(verts[i].xyz[1] > maxs.y)
           maxs.y = verts[i].xyz[1];

        lodbase.vects[i].xyz[2] = verts[i].xyz[2];
        if(verts[i].xyz[2] < mins.z)
           mins.z = verts[i].xyz[2];
        if(verts[i].xyz[2] > maxs.z)
           maxs.z = verts[i].xyz[2];

        lodbase.vects[i].normal[0] = verts[i].normal[0];
        lodbase.vects[i].normal[1] = verts[i].normal[1];
        lodbase.vects[i].normal[2] = verts[i].normal[2];
    }

    lodbase.mats.resize(materials_count);
    for(i=0; i<materials_count; i++)
    {
        lodbase.mats[i].diffuse[0] = materials[i].diffuse[0];
        lodbase.mats[i].diffuse[1] = materials[i].diffuse[1];
        lodbase.mats[i].diffuse[2] = materials[i].diffuse[2];
        lodbase.mats[i].specular[0] = materials[i].specular[0];
        lodbase.mats[i].specular[1] = materials[i].specular[1];
        lodbase.mats[i].specular[2] = materials[i].specular[2];
    }

    lodbase.tris.resize(trcount);
    for(i=0; i<trcount; i++)
    {
        lodbase.tris[i].verts[0] = triangles[i].verts[0];
        lodbase.vects[lodbase.tris[i].verts[0]].triangles.push_back(i);
        lodbase.tris[i].verts[1] = triangles[i].verts[1];
        lodbase.vects[lodbase.tris[i].verts[1]].triangles.push_back(i);
        lodbase.tris[i].verts[2] = triangles[i].verts[2];
        lodbase.vects[lodbase.tris[i].verts[2]].triangles.push_back(i);
    }

    lodbase.vcount = verts_count;
    lodbase.mcount = materials_count;
    lodbase.tcount = trcount;

    for(i=0; i<verts_count; i++)
       lodbase.vects[i].tcount = lodbase.vects[i].triangles.size();

    return 0;
}




float surface_factor(std::vector<tvertex>& verts, int vcount,
             std::vector<triangle>& triangles, int trcount,
             std::map<int,float>& trsurfaces,
             int v1, int v2)
{
    float actsurface = 0.f;
    float futuresurf = 0.f;
    int i;
    int end;
    vertex future_vert;
    triangle* tr;
    int v1idx, v2idx;
    int vcontained;
    std::map<int,float>::const_iterator it;
    float surf;

    average_vertices(verts[v1], verts[v2], future_vert);

    end = verts[v1].triangles.size();
    for(i=0; i<end; i++)
    {
        tr = &triangles[verts[v1].triangles[i]];
        vcontained = 0;

        if(tr->verts[0]==v1) {
            vcontained++;
            v1idx = 1;
            v2idx = 2;
        }
        else if(tr->verts[0]==v2) {
            vcontained++;
            v1idx = 2;
            v2idx = 1;
        }

        if(tr->verts[1]==v1) {
            vcontained++;
            v1idx = 0;
            v2idx = 2;
        }
        else if(tr->verts[1]==v2) {
            vcontained++;
            v1idx = 2;
            v2idx = 0;
        }

        if(tr->verts[2]==v1) {
            vcontained++;
            v1idx = 0;
            v2idx = 1;
        }
        else if(tr->verts[2]==v2) {
            vcontained++;
            v1idx = 1;
            v2idx = 0;
        }


        // Thanks to keeping already calculated triangles' surfaces,
        // we are not calculating them again.
        it = trsurfaces.find(verts[v1].triangles[i]);
        if(it!=trsurfaces.end())
            actsurface += (*it).second;
        else
        {
            surf = triangle_surface(verts, tr->verts[0], tr->verts[1], tr->verts[2]);
            trsurfaces[verts[v1].triangles[i]] = surf;
            actsurface += surf;
        }
        

        if(vcontained==1)
        {
            futuresurf += triangle_surface(verts, tr->verts[v1idx], tr->verts[v2idx], future_vert);
        }
        // else if(vcontained==2) // the triangle is being deleted
    }

    end = verts[v2].triangles.size();
    for(i=0; i<end; i++)
    {
        tr = &triangles[verts[v2].triangles[i]];
        vcontained = 0;

        if(tr->verts[0]==v1) {
            vcontained++;
            v1idx = 1;
            v2idx = 2;
        }
        else if(tr->verts[0]==v2) {
            vcontained++;
            v1idx = 2;
            v2idx = 1;
        }

        if(tr->verts[1]==v1) {
            vcontained++;
            v1idx = 0;
            v2idx = 2;
        }
        else if(tr->verts[1]==v2) {
            vcontained++;
            v1idx = 2;
            v2idx = 0;
        }

        if(tr->verts[2]==v1) {
            vcontained++;
            v1idx = 0;
            v2idx = 1;
        }
        else if(tr->verts[2]==v2) {
            vcontained++;
            v1idx = 1;
            v2idx = 0;
        }


        // This time we are not adding surfaces of triangles that contain both 
        // vertices (they have already been added).
        
        if(vcontained==1)
        {
            it = trsurfaces.find(verts[v2].triangles[i]);
            if(it!=trsurfaces.end())
                actsurface += (*it).second;
            else
            {
                surf = triangle_surface(verts, tr->verts[0], tr->verts[1], tr->verts[2]);
                trsurfaces[verts[v2].triangles[i]] = surf;
                actsurface += surf;
            }
            futuresurf += triangle_surface(verts, tr->verts[v1idx], tr->verts[v2idx], future_vert);
        }
        // else if(vcontained==2) // the triangle is being deleted
    }

    return futuresurf-actsurface;
}

// Based on a mesh surface change.
// Choice an edge that after collapsing would least change the surface of the mesh.
std::pair<int,int> find_edge_to_collapse(std::vector<tvertex>& verts, int vcount,
        std::vector<triangle>& triangles, int trcount)
{
    // Indicate, whether an edge has been already checked.
    std::map<int,bool> checked_edges;
    int i,j;
    int end;
    triangle* tr;
    float surface_mistake;
    float min_surface_mistake = 10e9;
    std::pair<int,int> min_edge(0,0);
    std::map<int,bool>::const_iterator endi;
    std::map<int,bool>::const_iterator it;
    std::list<float> mistakes;
    std::map<int,float> trsurfaces;
    for(i=0; i<vcount; i++)
    {
        // This vertex has been deleted.
        end = verts[i].triangles.size();
        if(end == 0)
            continue;

        // There are two possible problems decreasing speed of the algorithm:
        // - every edge has two endings (double checking of every edge)
        // - different triangles might consist the same edge (multiple checking of the same edge)
        // As a solution for the first, I'm always checking an edge (a,b), when a < b.
        // For the second one, I'm checking the edge only if it hasn't been checked already.
        checked_edges.clear();
        for(j=0; j<end; j++)
        {
            tr = &triangles[verts[i].triangles[j]];
            if(i < tr->verts[0])
            {
                endi = checked_edges.end();
                it = checked_edges.find(tr->verts[0]);
                if(it==endi)
                {
                    surface_mistake = surface_factor(verts, vcount, triangles, trcount, 
                              trsurfaces, i, tr->verts[0]);
                    mistakes.push_back(surface_mistake);
                    if(std::abs(surface_mistake) < std::abs(min_surface_mistake)) {
                        min_surface_mistake = surface_mistake;
                        min_edge.first = i;
                        min_edge.second = tr->verts[0];
                    }
                    checked_edges[tr->verts[0]] = true;
                }
            }
            if(i < tr->verts[1])
            {
                endi = checked_edges.end();
                it = checked_edges.find(tr->verts[1]);
                if(it==endi)
                {
                    surface_mistake = surface_factor(verts, vcount, triangles, trcount, 
                              trsurfaces, i, tr->verts[1]);
                    mistakes.push_back(surface_mistake);
                    if(std::abs(surface_mistake) < std::abs(min_surface_mistake)) {
                        min_surface_mistake = surface_mistake;
                        min_edge.first = i;
                        min_edge.second = tr->verts[1];
                    }
                    checked_edges[tr->verts[1]] = true;
                }
            }
            if(i < tr->verts[2])
            {
                endi = checked_edges.end();
                it = checked_edges.find(tr->verts[2]);
                if(it==endi)
                {
                    surface_mistake = surface_factor(verts, vcount, triangles, trcount, 
                              trsurfaces, i, tr->verts[2]);
                    mistakes.push_back(surface_mistake);
                    if(std::abs(surface_mistake) < std::abs(min_surface_mistake)) {
                        min_surface_mistake = surface_mistake;
                        min_edge.first = i;
                        min_edge.second = tr->verts[2];
                    }
                    checked_edges[tr->verts[2]] = true;
                }
            }
        }
    }
    //std::list<float>::iterator iter = mistakes.begin();
    //for(; iter!=mistakes.end(); iter++)
     //   cout << *iter << endl;
    cout << "minimum mistake = " << min_surface_mistake << std::endl;
    return min_edge;
}


std::pair<int,int> random_edge(std::vector<triangle>& triangles, int trcount)
{
    struct timeval t;
    int ret;
    int tr;
    int vert;
    int v[2];
    int i;
    // Initialize random number generator.
    ret = gettimeofday(&t, NULL);
    assert(ret==0);
    srandom((unsigned int)(t.tv_sec*1e6+t.tv_usec));
    tr = random()%trcount;
    vert = random()%3; // 3 vertices make a triangle
    i=0;
    if(vert!=0)
        v[i++] = triangles[tr].verts[0];
    if(vert!=1)
        v[i++] = triangles[tr].verts[1];
    if(vert!=2)
        v[i++] = triangles[tr].verts[2];
    assert(i==2);
    
    if(v[0]<v[1])
        return std::pair<int,int>(v[0],v[1]);
    else
        return std::pair<int,int>(v[1],v[0]);
}



// todo:
// optimize half-edge collapse (now it works same as full-collapse)
// full collapsing is based on the average between two vertices
// it's not checked if the newly created vertex by chance already exist
int edge_collapse(std::vector<tvertex>& verts, int& vcount, 
                std::vector<triangle>& triangles, int& trcount, int v1, int v2,
                int cpsetype)//collapsion type: half=0 or full=1
                //std::vector<int>& changedtri, std::vector<int>& delv)
                          // triangles changed     and vertices deleted
                          // during collapsion
{
    int i,end;
    triangle* tr;
    int v;
    int vidx;
    int vcontained;
    tvertex vert;
    int deltrcount = 0;
    int delvcount = 0;
    int v2idx;

    
    // Calculate the lower index.
    v = v1<v2?v1:v2;
    if(cpsetype==1)//full
       average_vertices(verts[v1], verts[v2], verts[v]);//put the new v at the lower idx
    //in half-collapsing I'm picking the vertex at the lower index

    end = verts[v1].triangles.size();
    for(i=0; i<end; i++)
    {
        // Find out if the triangle contains one of the vs or both.
        tr = &triangles[verts[v1].triangles[i]];
        vcontained = 0;
        if(tr->verts[0]==v1) {
            vidx=0;
            vcontained++;
        }
        else if(tr->verts[0]==v2) {
            v2idx = 0;
            vcontained++;
        }

        if(tr->verts[1]==v1) {
            vidx=1;
            vcontained++;
        }
        else if(tr->verts[1]==v2) {
            v2idx = 1;
            vcontained++;
        }

        if(tr->verts[2]==v1) {
            vidx=2;
            vcontained++;
        }
        else if(tr->verts[2]==v2) {
            v2idx = 2;
            vcontained++;
        }

        // Since triangles containing collapsing edge are deleted,
        // we leave them in peace (the vector of triangles is no longer coherent).
        if(vcontained==1)
        {
            // Update the list of vs of the t
            tr->verts[vidx] = v;
            // Memorize that the v belongs to the t
            vert.triangles.push_back(verts[v1].triangles[i]);
        }
        else// Both vs belong to the t
        {
           // Recognise which v is the 3rd one from the t
            if(vidx!=0 && v2idx!=0)
                vidx = 0;
            else if(vidx!=1 && v2idx!=1)
                vidx = 1;
            else if(vidx!=2 && v2idx!=2)
                vidx = 2;

            vidx = tr->verts[vidx];
            std::vector<int>::iterator j, endj;
            endj = verts[vidx].triangles.end();
            for(j=verts[vidx].triangles.begin(); j!=endj; j++)
                if((*j)==verts[v1].triangles[i])
                    break;
            assert(j!=endj);
            // Delete the t from the v's list of triangles
            verts[vidx].triangles.erase(j);
            verts[vidx].tcount--;
            if(verts[vidx].triangles.size()==0)
               // This was the last t of this v, we're are not keeping selfalone vs
                delvcount++;
            deltrcount++;
        }
    }

    end = verts[v2].triangles.size();
    for(i=0; i<end; i++)
    {
        tr = &triangles[verts[v2].triangles[i]];
        vcontained = 0;
        if(tr->verts[0]==v1)
            vcontained++;
        else if(tr->verts[0]==v2) {
            vidx=0;
            vcontained++;
        }

        if(tr->verts[1]==v1)
            vcontained++;
        else if(tr->verts[1]==v2) {
            vidx=1;
            vcontained++;
        }

        if(tr->verts[2]==v1)
            vcontained++;
        else if(tr->verts[2]==v2) {
            vidx=2;
            vcontained++;
        }

        if(vcontained==1)
        {
            tr->verts[vidx] = v;
            vert.triangles.push_back(verts[v2].triangles[i]);
        }
    }

    if(vert.triangles.size()>0)
        verts[v].triangles = vert.triangles;
    else {
        verts[v].triangles.clear();
        delvcount++;
    }
    

    if(v==v1)
        v=v2;
    else
        v=v1;
    // The v at the higher index is now free, mark it.
    verts[v].triangles.clear();
    delvcount++;

    trcount -= deltrcount;
    vcount -= delvcount;
    return 0;
}


int clear_obj_data(object3dv<tvertex>& dirty, object3dv<vertex>& newlod,
               int vcount, int trcount,//materials are copied
               vect& min, vect& max)
{
    std::map<int,int> move;
    int i, end;
    std::map<int,bool> validtris;
    std::map<int,bool>::const_iterator vtit;
    std::map<int,bool>::const_iterator vtendit;


    // Move the vertices to the correct places and save their positions
    // for future corrections of triangles' indices.
    newlod.vects.resize(vcount);
    i=0;
    end=dirty.vects.size();
    for(int j=0; j<end && i<vcount; j++)
        if(dirty.vects[j].triangles.size()!=0)
        {
            move[j] = i;
            newlod.vects[i]=dirty.vects[j];
            const vertex& v=newlod.vects[i];
            if(v.xyz[0]<min.x)
               min.x=v.xyz[0];
            if(v.xyz[0]>max.x)
               max.x=v.xyz[0];
            if(v.xyz[1]<min.y)
               min.y=v.xyz[1];
            if(v.xyz[1]>max.y)
               max.y=v.xyz[1];
            if(v.xyz[2]<min.z)
               min.z=v.xyz[2];
            if(v.xyz[2]>max.z)
               max.z=v.xyz[2];
            int k=0;
            int endk=dirty.vects[j].triangles.size();
            for(; k<endk; k++)
                validtris[dirty.vects[j].triangles[k]] = true;
            i++;
        }
    
    assert(i==vcount);
    assert((int)validtris.size()==trcount);

    newlod.tris.resize(trcount);

    vtendit=validtris.end();
    i=0;
    int j=0;
    for(vtit=validtris.begin(); vtit!=vtendit; vtit++) {
          // move the triangle to the correct place
          // correct vertices of this triangle
          newlod.tris[i].verts[0]=move[dirty.tris[(*vtit).first].verts[0]];
          newlod.tris[i].verts[1]=move[dirty.tris[(*vtit).first].verts[1]];
          newlod.tris[i].verts[2]=move[dirty.tris[(*vtit).first].verts[2]];
          
       i++;
       j++;
    }
    assert(i==trcount);

    newlod.vcount=vcount;
    newlod.tcount=trcount;
    newlod.mats=dirty.mats;
    newlod.mcount=dirty.mcount;
    
    return 0;
}
    

//static void disperse(int vdelnum, float xmin, float xmax, float ymin, float ymax,
      //float zmin, float zmax, int& xparts, int& yparts, int& zparts)
static void disperse(int vdelnum, const vect& mins, const vect& maxs, 
                     int& xparts, int& yparts, int& zparts)
{
   // calculate a distribution of variables x,y & z
   int parts[3];
   float tmp;
   
   tmp = maxs.x-mins.x + maxs.y-mins.y + maxs.z-mins.z;
   parts[0] = (int)std::pow(vdelnum,(maxs.x-mins.x)/tmp);
   parts[1] = (int)std::pow(vdelnum,(maxs.y-mins.y)/tmp);
   parts[2] = (int)std::pow(vdelnum,(maxs.z-mins.z)/tmp);

   // extend the number of grids as less as possible
   if(parts[0]*parts[1]*parts[2] < vdelnum) {
      bool posused[3] = {false,false,false};
      int idx;
      if(parts[0]<parts[1]) {
         tmp = parts[0];
         idx = 0;
      }
      else {
         tmp = parts[1];
         idx = 1;
      }
      
      if(parts[2]<tmp) {
         tmp = parts[2];
         idx = 2;
      }
      
      parts[idx]++;
      posused[idx]=true;

      if(parts[0]*parts[1]*parts[2] < vdelnum) {
         int idx2;

         if(posused[0]) {
            idx = 1;
            idx2 = 2;
         }
         else if(posused[1]) {
            idx = 0;
            idx2 = 2;
         }
         else {
            idx = 0;
            idx2 = 1;
         }

         if(parts[idx]<parts[idx2]) {
            parts[idx]++;
            posused[idx]=true;
         }
         else {
            parts[idx2]++;
            posused[idx2]=true;
         }

         // todo
         // will it ever happen?
         if(parts[0]*parts[1]*parts[2] < vdelnum) {
            if(posused[0]==false)
               parts[0]++;
            else if(posused[1]==false)
               parts[1]++;
            else
               parts[2]++;
         }
      }
   }

   xparts = parts[0];
   yparts = parts[1];
   zparts = parts[2];
}

// Make LOD from level <level> to level <tolevel>.
// The LOD on level <level> should be put to lodbases[tolevel], 
// the function will change it making new LOD.
// todo
// a copy of all of the vertices are being made 
// to calculate the split points defining the grids
int dispersed_cubes_lod(std::vector< object3dv<tvertex> >& lodbases, 
         const std::vector< std::pair<vect,vect> >& bases_ranges,
         int level, int tolevel, int vdelnum, int edgecolltype,
         int& newvcount, int& newtcount)
{
   int xparts, yparts, zparts;
   std::vector<float> v;
   std::vector<float> splits[3];
   int i;
   int end;
   // grid structure 
   std::vector<int>** grid;

   disperse(vdelnum, (bases_ranges[level]).first, (bases_ranges[level]).second,
            xparts,yparts,zparts);
   
   object3dv<tvertex>& lodb = lodbases[tolevel];
   std::vector<tvertex>& vects = lodb.vects;

   // the split points calculated below are being made unprecisely
   // to be precise we have to find y middle points seperately for each x-slice
   // and then again inside each xy-dice
   // x
   end = lodb.vcount;
   v.resize(lodb.vcount);
   for(i=0; i<end; i++)
      v[i] = vects[i].xyz[0];
   std::sort(v.begin(),v.end());

   // splits define the shape of the grid
   splits[0].resize(xparts-1);
   end = xparts-1;
   for(i=0; i<end; i++) {
      //cout<<((i+1)*lodb.vcount/xparts)<<endl;
      splits[0][i] = v[(i+1)*lodb.vcount/xparts];
   }
   

   // y
   for(i=0; i<end; i++)
      v[i] = vects[i].xyz[1];
   std::sort(v.begin(),v.end());

   splits[1].resize(yparts-1);
   end = yparts-1;
   for(i=0; i<end; i++) {
      //cout<<((i+1)*lodb.vcount/yparts)<<endl;
      splits[1][i] = v[(i+1)*lodb.vcount/yparts];
   }

   // z
   for(i=0; i<end; i++)
      v[i] = vects[i].xyz[2];
   std::sort(v.begin(),v.end());

   splits[2].resize(zparts-1);
   end = zparts-1;
   for(i=0; i<end; i++) {
      //cout<<((i+1)*lodb.vcount/zparts)<<endl;
      splits[2][i] = v[(i+1)*lodb.vcount/zparts];
   }

   // allocating a table of vectors makes troubles
   // it seems that destructors of the vectors are being called
   // after deallocation of the vectors' table
   grid = new std::vector<int>*[xparts*yparts*zparts];
   end=xparts*yparts*zparts;
   for(i=0; i<end; i++)
      grid[i] = new std::vector<int>();

   end = lodb.vcount;
   int j,endj;
   int k,endk;
   int l,endl;
   endj = xparts-1;
   for(i=0; i<end; i++) {
      for(j=0; j<endj; j++) {
         if(lodb.vects[i].xyz[0] <= splits[0][j])
            break;
      }
      endk=yparts-1;
      for(k=0; k<endk; k++) {
         if(lodb.vects[i].xyz[1] <= splits[1][k])
            break;
      }
      endl=zparts-1;
      for(l=0; l<endl; l++) {
         if(lodb.vects[i].xyz[2] <= splits[2][l])
            break;
      }
           //j*yparts*zparts + k*zparts + l
      (*grid[(j*yparts+k)*zparts + l]).push_back(i);
   }
   
   /*for(j=0; j<xparts; j++) {
      if(j<xparts-1)
         cout<<splits[0][j]<<endl;
      for(k=0; k<yparts; k++) {
         for(l=0; l<zparts; l++) {
            std::vector<int>& ve = (*grid[(j*yparts+k)*zparts + l]);
            end=ve.size();
            for(i=0; i<end; i++)
               cout<<(lodb.vects[ve[i]].xyz[0])<<endl;
         }
      }
   }*/
   
   // finally, lets proceed with some lodding...
   j=0;//grid iterator
   int mintcount;//minimal number of triangles containing vertex l
   l=0;//vertex l
   int gridcount=xparts*yparts*zparts;
   int lastvisited=-1;
   int whiledog=-1;
   int vcount=lodb.vcount;
   int tcount=lodb.tcount;
   std::map<int,bool> delv;
   std::map<int,bool>::iterator it;
   
   for(i=0; i<vdelnum;)
   {
      mintcount=(int)INF;
      while((*grid[j]).size()<=0)
      {
         if(whiledog==-1)
            whiledog=j;
         else if(whiledog==j) {
            break;
         }
         j++;
         if(j==gridcount)
            j=0;
      }
      if(whiledog==j) {
         //todo doesn't work for sphere.gat
         //std::cerr<<"break1"<<std::endl;
         break;
      }
      else
         whiledog=-1;

      end=(*grid[j]).size();
      for(k=0; k<end; k++)
      {
         // todo
         // it might be done smarter, 
         // e.g. assuming that the center of the object 
         // lies in the center of the coordinate system, 
         // we can start from the outest grids 
         // moving to the inside of the system 
         if(lodb.vects[(*grid[j])[k]].tcount < mintcount) {
            mintcount=lodb.vects[(*grid[j])[k]].tcount;
            l=k;
         }
      }
      if(mintcount >= TRIANGLES_CONSTRAINT) {
         if(lastvisited!=j) {
            // go to the next grid
            if(lastvisited==-1)
               lastvisited=j;//mark the first grid that we jumped over 
                             //because of the TRIANGLES_CONSTRAINT
            j++;
            if(j>=gridcount)
               j=0;
            //i--;
            continue;
         }
         else {//we made a circle and came back to the same grid
            lastvisited=-1;
         }
      }
      
      // we got a grid with a satisfactory vertex
      // try to pick up the other end of the edge from the same grid
      bool pickif=true;
      int l2;
      if((*grid[j]).size()>1) {
         for(k=0; k<end; k++)
         {
            if(k!=l && in_one_triangle(lodb.vects, lodb.tris,
                                       (*grid[j])[l],(*grid[j])[k])) {
               if(lodb.vects[(*grid[j])[k]].tcount > TRIANGLES_CONSTRAINT)
                  continue;
               
               // don't search for minimal connections with the second vertex
               pickif=false;
               l2=(*grid[j])[k];
               break;
            }
         }
      }
      //todo
      if(pickif)
      {
         // pick the second v of the edge from the first common triangle
         // todo
         // minimal connections criteria can be also applied here
         // moreover, it can be combined with the surface change criteria
         // e.g. a*surface_change + b*minimal_connections
         const std::vector<int>& t1 = lodb.vects[(*grid[j])[l]].triangles;
         /*if(lodb.tris[t1[0]].verts[0] != (*grid[j])[l])
            l2 = lodb.tris[t1[0]].verts[0];
         else
            l2 = lodb.tris[t1[0]].verts[1];*/
         endj=lodb.vects[(*grid[j])[l]].tcount;
         int m1,m2;
         int jj;
         for(jj=0; jj<endj; jj++)
         {
            // find an index of v l
            if(lodb.tris[t1[jj]].verts[0] == (*grid[j])[l]) {
               m1=lodb.tris[t1[jj]].verts[1];
               m2=lodb.tris[t1[jj]].verts[2];
            }
            else if(lodb.tris[t1[jj]].verts[1] == (*grid[j])[l]) {
               m1=lodb.tris[t1[jj]].verts[0];
               m2=lodb.tris[t1[jj]].verts[2];
            }
            else {
               m1=lodb.tris[t1[jj]].verts[0];
               m2=lodb.tris[t1[jj]].verts[1];
            }

            it=delv.find(m1);
            if(it==delv.end()) {//the neighbourhood of the v m1 has not been changed
               l2=m1;
               break;
            }
            it=delv.find(m2);
            if(it==delv.end()) {
               l2=m2;
               break;
            }
         }

         //assert(jj<endj);
         if(jj>=endj)
         {           //all of the neighbour vertices of v l
                     //were affected in former lodding, so this v
                     //belongs to no edge in the grid structure;
                     //delete it
            int kk,ll,endkk,n;
            std::vector<int>::iterator m;
            std::vector<int>::iterator endm;
            int lv=(*grid[j])[l];

            endj = xparts-1;
            for(jj=0; jj<endj; jj++) {
               if(lodb.vects[lv].xyz[0] <= splits[0][jj])
                  break;
            }
            endkk = yparts-1;
            for(kk=0; kk<endkk; kk++) {
               if(lodb.vects[lv].xyz[1] <= splits[1][kk])
                  break;
            }
            endl = zparts-1;
            for(ll=0; ll<endl; ll++) {
               if(lodb.vects[lv].xyz[2] <= splits[2][ll])
                  break;
            }
            n=(jj*yparts+kk)*zparts + ll;
            endm=(*grid[n]).end();
            for(m=(*grid[n]).begin(); m<endm; m++) {
               if((*m) == lv) {
                  (*grid[n]).erase(m);
                  break;
               }
            }
            assert(m<endm);
            delv[lv]=true;
            //i--;
            continue;
         }
      }

      l=(*grid[j])[l];
      //l2=(*grid[j])[l2];
      
      
      // Since this algorithm is supposed to disperse changes,
      // to collapse I'm picking only edges from isolated regions
      // (this way collapsions do not affect each other and they are also scattered).
      {
         int jj,kk,ll,endkk;

         int n;
         std::vector<int>::iterator m;
         std::vector<int>::iterator endm;
         // delete v l from the grid structure (don't pick it up again)
         endj = xparts-1;
         for(jj=0; jj<endj; jj++) {
            if(lodb.vects[l].xyz[0] <= splits[0][jj])
               break;
         }
         endkk = yparts-1;
         for(kk=0; kk<endkk; kk++) {
            if(lodb.vects[l].xyz[1] <= splits[1][kk])
               break;
         }
         endl = zparts-1;
         for(ll=0; ll<endl; ll++) {
            if(lodb.vects[l].xyz[2] <= splits[2][ll])
               break;
         }
         n=(jj*yparts+kk)*zparts + ll;
         endm=(*grid[n]).end();
         for(m=(*grid[n]).begin(); m<endm; m++) {
            if((*m) == l) {
               (*grid[n]).erase(m);
               break;
            }
         }
         assert(m<endm);
         delv[l]=true;//it's not necessary, just to be sure


         // delete v l2 from the grid structure
         for(jj=0; jj<endj; jj++) {
            if(lodb.vects[l2].xyz[0] <= splits[0][jj])
               break;
         }
         for(kk=0; kk<endkk; kk++) {
            if(lodb.vects[l2].xyz[1] <= splits[1][kk])
               break;
         }
         for(ll=0; ll<endl; ll++) {
            if(lodb.vects[l2].xyz[2] <= splits[2][ll])
               break;
         }
         n=(jj*yparts+kk)*zparts + ll;
         endm=(*grid[n]).end();
         for(m=(*grid[n]).begin(); m<endm; m++)
         {
            if((*m) == l2) {
               (*grid[n]).erase(m);
               break;
            }
         }
         assert(m<endm);
         delv[l2]=true;//it's not necessary, just to be sure

         
         int ve;
         triangle* trp;
         endk=lodb.vects[l].tcount;
      
         // delete neighbourhood of the v l
         for(k=0; k<endk; k++) {
            trp = &lodb.tris[lodb.vects[l].triangles[k]];
            //todo
            //it's not worth to make such a small loop
            for(int w=0; w<3; w++)
            {
               ve = (*trp).verts[w];
               if(ve!=l && ve!=l2) {
                  it=delv.find(ve);
                  if(it!=delv.end())
                     continue;//v has been already deleted

                  // now comes a tricky step: find the grid to which belongs vertex ve
                  for(jj=0; jj<endj; jj++) {
                     if(lodb.vects[ve].xyz[0] <= splits[0][jj])
                        break;
                  }
                  for(kk=0; kk<endkk; kk++) {
                     if(lodb.vects[ve].xyz[1] <= splits[1][kk])
                        break;
                  }
                  for(ll=0; ll<endl; ll++) {
                     if(lodb.vects[ve].xyz[2] <= splits[2][ll])
                        break;
                  }
                  
                  n=(jj*yparts+kk)*zparts + ll;
                  endm=(*grid[n]).end();
                  for(m=(*grid[n]).begin(); m<endm; m++) {
                     if((*m) == ve) {
                        (*grid[n]).erase(m);
                        break;
                     }
                  }
                  assert(m<endm);

                  delv[ve]=true;
               }
            }
            
         }

         // delete neighbourhood of the v l2
         endk=lodb.vects[l2].tcount;
         for(k=0; k<endk; k++) {
            trp = &lodb.tris[lodb.vects[l2].triangles[k]];
            //todo
            //it's not worth to make such a small loop
            for(int w=0; w<3; w++) {
               ve = (*trp).verts[w];
               if(ve!=l && ve!=l2) {
                  it=delv.find(ve);
                  if(it!=delv.end())
                     continue;//v has been already deleted

                  // now comes a tricky step: find the grid to which belongs vertex ve
                  for(jj=0; jj<endj; jj++) {
                     if(lodb.vects[ve].xyz[0] <= splits[0][jj])
                        break;
                  }
                  for(kk=0; kk<endkk; kk++) {
                     if(lodb.vects[ve].xyz[1] <= splits[1][kk])
                        break;
                  }
                  for(ll=0; ll<endl; ll++) {
                     if(lodb.vects[ve].xyz[2] <= splits[2][ll])
                        break;
                  }
                  
                  n=(jj*yparts+kk)*zparts + ll;
                  endm=(*grid[n]).end();
                  for(m=(*grid[n]).begin(); m<endm; m++) {
                     if((*m) == ve) {
                        (*grid[n]).erase(m);
                        break;
                     }
                  }
                  assert(m<endm);

                  delv[ve]=true;
               }
            }
            
         }
      }
      
      int n=tcount;
      // uff..., after all, collapse the edge (l,l2)
      edge_collapse(lodb.vects, vcount, lodb.tris, tcount,
                  l, l2, edgecolltype);

      //std::cerr<<tcount<<std::endl;
      i+=n-tcount;
      
   }
   cerr<<"to"<<vdelnum<<std::endl;
   cerr<<"aa"<<i<<std::endl;
   //cout<<"v "<<lodb.vects.size()<<std::endl;
   //cout<<"t "<<lodb.tris.size()<<std::endl;
   //clear_obj_data(lodb.vects,vcount,lodb.tris,tcount);
   newvcount=vcount;
   newtcount=tcount;
   //lodb.vcount=vcount;
   //lodb.tcount=tcount;
   //cout<<"v "<<lodb.vects.size()<<std::endl;
   //cout<<"t "<<lodb.tris.size()<<std::endl;
   

   // deallocation of the grid structure
   end=xparts*yparts*zparts;
   for(i=0; i<end; i++)
      delete grid[i];
   delete grid;

   return 0;
}




// Helper function for distant_edges()
void recursive_distant_edges(std::vector<tvertex>& verts, int vcount, 
        std::vector<int>& result, float& mindistance, int my_num, int istart, 
        std::vector<int>& chosenv)
{
    int i;
    if(my_num==0)
    {
        for(i=istart; i<vcount; i++)
        {
            //for(int j=chosenv.size()-1; j>0; j--)
             //   cout<<chosenv[j]<<' ';
            //cout<<i<<endl;
            mindistance++;
        }
    }
    // make recursive call
    else if(istart+1<vcount)
    {
        for(i=istart; i<vcount; i++)
        {
            chosenv[my_num] = i;
            recursive_distant_edges(verts,vcount,result,mindistance,my_num-1,i+1,chosenv);
        }
    }
}

/*vector<std::pair<int,int> >*/void distant_edges(std::vector<tvertex>& verts, int& vcount, 
                std::vector<triangle>& triangles, int& trcount, int n)
{
    assert(n>0);
    //bool closer = false;
    std::vector<int> result;
    float mindistance;
    std::vector<int> chosenv;
    
    chosenv.resize(n);
    result.resize(n);
    mindistance=1e12;
    mindistance=0;
    recursive_distant_edges(verts, vcount, result, mindistance, n-1, 0,chosenv);
    cout<<mindistance<<std::endl;
}


bool in_one_triangle(const std::vector<tvertex>& verts,
         const std::vector<triangle>& tris, int v1, int v2)
{
   const std::vector<int>& t1 = verts[v1].triangles;
   int end=verts[v1].tcount;
   for(int i=0; i<end; i++) {
      if(tris[t1[i]].verts[0] == v2)
         return true;
      else if(tris[t1[i]].verts[1] == v2)
         return true;
      else if(tris[t1[i]].verts[2] == v2)
         return true;
   }
   return false;
}
